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Abstract 

We investigate the hetero-epitaxial growth of thin film deposited on a (001) substrate via molec- 
ular dynamics simulations, using six fee transition metals as our modeling systems. By studying 
the radial distribution function in the film layers, we demonstrate the importance of the sign of 
lattice mismatch on the layer structure. For positive lattice mismatches, the film favors pseudo- 
morphic growth, whereas for negative mismatches, a sharp transition happens within the first few 
monolayers of adatoms and the film layers are transformed into the close-packed (111) structure. 
We propose a way to quantify the compositional percentage of different planar structures in an 
epilayer, and demonstrate the strong asymmetric effect between the tensile and compressive cases 
of deposition. How temperature affects the asymmetry is also discussed. 

PACS numbers: 68.55.-a, 81.15.Aa, 68.47.De, 68.35.Ct 
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How to epitaxially grow thin films on top of a substrate of different crystalline structure 
is a challenging task and continually attracts great attention of researchers due to its wide 
applications^. It is commonly acknowledged that to grow a high-quality single-crystal film 
in a system where the lattice constant of the grown material and the substrate differs is 
very difficult. Therefore, a small lattice mismatch is generally recognized as a determina- 
tive factor for epitaxial growth, together with other factors such as surface/interface free 
energy, chemical bondings, temperature, etc.—. In hetero-epit axial growth, an elastic strain 
is built up due to mismatched structure. Misfit dislocation is thus formed in the film lay- 
ers to release the elastic energy while the elastic strain approaches a certain degree, which 
allows us to define a critical thickness of epitaxy. In the classical continuum theory^^, 
the elastic energy is simply determined by the square of the lattice mismatch defined by 
/ = (a s — af)/a,f where a/ and a s are the unstrained lattice constants of the film and the 
substrate, respectively. As a consequence, the critical thickness depend only on |/|. This is 
why the epitaxial growth of thin film was usually examined regardless of the sign of f—^-. 
Nevertheless, scientists have shown that the sign of mismatch does have a significant effect 
on the structure and the formation of misfit dislocations, demonstrated by experiments and 
numerical studies*^ 1 ^. This effect can be attributed to the tensile-compressive asymmetry, 
originated from the anharmonicity of the atomic interaction^. Furthermore, surface energy 
also plays an important role in the film growth and roughening. It is known that the (111) 
surface has the lowest free energy for face-centered cubic (fee) materials. Therefore, an 
epitaxial growth on the (100) surface for fee metals will compete with the formation of the 
(111) structure. Different kinds of film structure growing on a (100) substrate have been in- 
vestigated in literatures^ 1 ^ 1 ^. Nonetheless, a systematic study of the positive and negative 
lattice mismatch on the structure of a film layer has rarely been conducted. The information 
is missed: to what degree can the tensile-compressive anharmonicity break the symmetry of 
the epitaxial growth predicted by the continuum theory? To answer this question, we chose 
six transition metals of fee crystalline structure as our modeling systems and studied epitax- 
ial growth of one metal on the other by means of molecular dynamics (MD) simulations. By 
using the radial distribution function, we were able to calculate the percentage of different 
surface structure in a film layer and studied the effect of lattice mismatch quantitatively. 

The six transition metals that we chose are the elements Ni(3.52A), Pd(3.89A), Pt(3.92A), 
in the group 8B, and Cu(3.6lA), Ag(4.09A), Au(4.08A), in the group IB. The corresponding 
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lattice constant^ has been given, following the element symbol. We studied the thin film 
deposition of a metal F onto a metal S where F and S were chosen from the six elements. 
There are totally 36 combinations, which give the value of lattice mismatch / ranging from 
— 14% to 16%. We studied the film growth on the (001) plane of a substrate and denoted 
the system by the notation "F/S(001)". The simulations were performed using a modified 
MD package, LAMMPS 17 , with a time step equal to 0.001 picosecond (ps). The interaction 
between atoms is modeled by the embedded-atom method 18 . The size of the substrate is 
of 10 unit cells in the [100] (x-) and [010] (y-) directions; thus, 200 atoms constitutes an 
atomic layer of the substrate. There are 9 atomic layers in a substrate. The atoms in the 
two bottom layers are fixed on the lattice positions. Above them, six layers are thermal 
layers subject to a temperature control by velocity rescaling method. The top layer of the 
substrate is free of the temperature control. Four temperatures were studied: 80K, 300K, 
600K and 900K. In order to simulate an infinite surface, periodic boundary condition was 
applied in x- and ^-directions. 

In simulating the depositional growth, the adatoms were randomly dropped, one by one, 
on the top of the substrate surface, 100A height above which, with an initial downward 
velocity corresponding to an incident energy 0.1 eV. This incident energy stays in the range 
of a typical deposition by evaporation. The deposition rate is one atom per 5 ps. Since 
one monolayer (ML) contains Nml = 200 atoms in our study, the growth rate is 1 ML 
per nanosecond, which is several orders of magnitude higher than in molecular beam epi- 
taxy. This choice is due to the limitation of standard MD techniques, restricted by today's 
computing power. However, the high deposition rate used here does not invalidate the 
present study. In laser deposition techniques, such as those used in Ref.— , the fluxes are 
relatively high, approaching to our setup. To clarify the effect of high deposition rate, we 
have measured the relaxation time for an adatom after collision with the surface. We found 
that it is less than 1 ps, consistent with theoretical predictions 19 . Therefore, the atom has 
enough time to relax to, at least, a local free energy minimum before coming in the next 
atom by this deposition rate. Another problem with high deposition rate is that there is no 
time for activation of many interlayer processes in film growth, such as exchange diffusion 
and multiple-atom concerted diffusion 20 . To take into account these infrequent complex 
processes, more sophisticated method such as hyperdynamics or kinetic Monte Carlo simu- 
lations should be applied 21 , where the flux can be set to a few orders of magnitude smaller 
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than in standard MD simulations. Nevertheless, it requires much more computing resources 
and efforts. For a systematic study of the asymmetric effect in heteroepitaxial growth, we 
neglected these infrequent processes as the first-order approximation and approached this 
problem using standard MD simulations, in order to capture the fundamental picture of such 
systems. The simulation results will be benchmarked later with experimental data obtained 
by other groups to verify the validness of this setup. We deposited 16 monolayers (ML) of 
adatoms for each studied case. In order to efficiently remove the heat brought in by the 
adatoms, the thermal layers were extended upward, one layer by one layer, each time when 
one ML of adatoms has been deposited onto the surface. 

We first studied the Cu/Ni system. The mismatch is —2.49%. The negative mismatch 
means that the film suffers from a compressive strain. We monitored the coverage of each 
epilayer during deposition by Ci(m) = Ne(m)/Nuh where Ne(m) is the number of adatoms 
in the £-th epilayer while m ML of adatom has been deposited. The results are shown in 
Fig. Q] at three temperatures: 80K, 300K and 600K. 




FIG. 1: Ci{m) and R(m) for Cu/Ni at 80K, 300K and 600K. The thin curves appearing from left 
to right at a given temperature denote, in turn, the coverage for £ = 1,2,... epilayer. The values 
are read from the left y-axis. R{m) is plotted in thick curve. The values are read from the right 
y-axis. 

Ci(m) is zero at beginning. It rapidly increases, one layer after another, and saturates 
to some value. For the first few layers, this value is nearly 100%, which indicates a good 
coverage. However, as £ increases, CV(m) saturates to a value significantly small then 1, 
suggesting the formation of voids in the high epilayers. We found that the higher the 
temperature, the less the voids will be formed. It can be attributed to the high kinetic 
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energy of the adatoms due to the high substrate temperature. With high kinetic energy, 
the adatoms can easily overcome energy barrier and migrate to the voids, flattening the film 
with good coverage. At 600K, the layer coverage is nearly 100% and no void is formed in 
the epilayers. Moreover, only a single Ce(m) curve increases at any given moment at this 
temperature. The film thus grows layer by layer; it is the Frank-van der Merwe growth mode. 
On the other hand, at 80K and 300K, following the layer-by-layer growth, several Ce(m) 
curves increases at the same time. The system undergoes a wetting-layer and island-growth 
mode, or Stranski-Krastanov growth mode. 

We also calculated the real-time roughness R(m) of the film surface. The definition of 
R(m) is the square root of the variance of the height of the atoms homogeneously sampling 
on the surface. The results are plotted in Fig. [T] in thick curves in unit of the thickness of 
a film layer. It clearly shows that the roughness is reduced by increasing temperature. At 
80K, R{m) is about 2.5^ after deposition of 16 ML, because of the island growth; but at 
600K, it is ca. 0.7d,L, owing to the flat-film growth. We have examined other systems, such 
as Pt/Cu (/ = —7.91%) and Pd/Ni (/ = —9.51%). The results showed similar trends of 
behavior: the coverage and the surface flatness are both improved by increasing the substrate 
temperature. 

We then investigated the structure of thin film after deposition by analyzing the two- 
dimensional radial distribution function (2D RDF) of each epilayer. The 2D RDF is calcu- 
lated by ge(r\\) = ae(r\\)/aeo, where <Ji{r\\) is the surface density of adatom in the £th epilayer 
at a radial distance r\\ away from an adatom, and = Ni/L 2 is the average surface density 
in the epilayer. In order to study the tensile-compressive asymmetry, ge(rn) for the two 
counter-systems, Pd/Cu (/ = —7.20%) and Cu/Pd (/ = 7.76%), at 300K were calculated. 
The results are illustrated in Fig. [2j 

Focus firstly on the tensile case, Cu/Pd. The referenced curve for the Pd substrate surface 
shows the standard 2D RDF for a (OOl)-plane in which the peaks appear subsequently at 
the positions a, V2a, 2a, ^/Ea, etc., with a being the nearest atomic distance. We have 
verified that the integration of these peaks with the appropriate weight 27rr||0"oo reproduces, 
respectively, the numbers of the nearest, the next nearest, the 3rd nearest neighbor atoms, 
and so forth. We found that the Cu epilayers display a gi(r\\) similar to the structure of Pd 
surface. It suggests a pseudomorphic growth of thin film along the [001]-direction. Since 
atomic spacing of Cu is smaller than Pd, the peaks in the Cu epilayers move slightly toward 
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FIG. 2: ^(rn) for Pd/Cu and Cu/Pd at 300K. The number on the left-hand side of a curve denotes 
the layer number I. For the reason of clarity, we have drawn the curves separately by shifting them 
vertically. ge(r\\) for the two substrate surfaces, Cu(001) and Pd(001), are plotted as references. 
The positions of peaks in <7io(j"n) are indicated, in turn, by a, \^2a, and so forth. 

left, as increasing the layer number, to reduce the strain. It approaches the bulk spacing 
2. 55 A of copper at the 10th epilayer. 

Focus secondly on the counter system, Pd/Cu. The Pd epilayers now suffer from com- 
pressive strain and release the strain by dislocation or alternating the film to a close-packed 
structure. Dramatic change in the ge(r\\) structure was observed, compared to the previous 
case. The change includes the change of the peak height and the appearance of some new 
peaks in reference of the 2D RDF of the Cu substrate. Note that a new peak appears at 
v3a since the 2nd epilayer. It strongly suggests a hexagonal arrangement of atoms because 
the typical 2D RDF for a hexagonal lattice shows a series of peaks at the positions a, \/3a, 
2a, \/7a, and so on. Therefore, the epilayers take a mixed structure combining the (001)- 
and the (lll)-planes. Recent experimental study 10 showed that the tensile overlayer for 
Cu/Pd(100) remains coherent up to about 9 ML, whereas for Pd/Cu(100), misfit disloca- 
tion appears just after few ML. To make a comparison with the experiments, we calculated 
in Fig.[3]the in-plane lattice strain (defined as (abuik - Ofiim)/abuik) and the interlayer distance 
of film during the growth process. 

We observed the striking asymmetric behavior in the relief of in-plane strain: the com- 
pressive strain relaxes within just few ML, much more quickly than the tensile strain does. 
This is in agreement with the experiments 10 . For interlayer distance, a sudden increase 
appears for the Pd/Cu system while the film grows more than one ML. This is because Pd 
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FIG. 3: (a) in-plane strain and (b) interlayer distance, during the film growth process 



has larger lattice constant than Cu; therefore, the distance between the Cu-Pd interface is 
smaller than the Pd-Pd epilayer distance inside the film. Similar argument can also explain 
the sudden decrease of the interlayer distance at one ML for the Cu/Pd system. Moreover, 
the interlayer distance inside the Pd film is about 2.lA. This value is larger than 1.945A, 
the layer distance of a bulk palladium along the [001] direction, but smaller than 2.246A, 
the distance along the [111] direction. Since the in-plane strain relaxes quickly in this case, 
this intermediate value reflects the fact of the formation of film layers, each of which takes a 
mixed planar (001)- and (lll)-structure. On the other hand, the interlayer distance inside 
the Cu film is smaller than 1.805A, half of the lattice constant of copper. Since Cu film grows 
pseudomorphically on Pd(001) substrate, this small value indicates the falling-down of an 
epilayer above another epilayer inside the film due to the unrelieved strain which has a in- 
plane lattice constant larger than the bulk one. The interlayer distance is therefore reduced. 
Scientists^ have observed exactly what we found here. Our results show good agreements 
with experiments in many aspects, more than qualitatively, although the deposition rate is 
fast. This gives a direct support to our simulations and the standard MD techniques can 
capture the main trend of behavior of this kind of systems and provide valuable information 
and insight into the problems. 

In order to have a direct image of the deposited layer structure in mind, we present in 
Fig. H] the snapshots of the 10th epilayer for the two systems. The Pd epilayer shows lattice 
displacement along the diagonal ([110]) direction to reduce compressive strain, resulting in 
a series of structural changes between (001)- and (lll)-planar structures, in form of strips. 
In the snapshot, some typical atomic distances, which contribute the peaks in the RDF 
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FIG. 4: Snapshots of the 10th epilayers for Pd/Cu and Cu/Pd at 300K. Some atomic distances 
which contribute the peaks in <%(rii) in Fig. [2] have been plotted and numbered. Dashed lines 
denote domain boundaries. 

in Fig. [21 are plotted and marked by the peak number. For example, the 2nd and the 
5th peaks in gt{r\\) are derived from the atoms in the (OOl)-planar strips, with the atomic 
arrangements similar to the links marked by 2 and 5, respectively, in the picture; the 3rd 
peak comes from the arrangement in the (lll)-planar strips, similar to the link 3. The 
links across the boundary of the (001)- and (lll)-planar strips contribute other peaks not 
appeared in the standard RDFs, such as the peak at rn = 2.4a which is due to the link 
6. The link 4' is across the strip boundary and has the length 1.93a. Nonetheless, lattice 
vibration at 300K smears the peak out, which becomes indistinguishable with the 4th peak 
at position 2a. 

The 10th epilayer for the Cu/Pd system shows much simpler structure, basically following 
the structure of the (OOl)-plane. It is an epitaxial growth. However, the epilayer feels tensile 
strain due to the Pd substrate. The strain is relaxed by creating breaks in the diagonal 
directions, separating the film layer into domains. This kind of rectangular domains has 
been experimentally observecr^ 1 ^. Moreover, we found that some atoms can sit on the 
breaks with the atom height deviating of the layer, which results in the small hump in 
the RDF curve near rii = 2.0A. For a link across a domain boundary, like the link 2', the 
length of the link is larger than the link inside a domain, like link 2. It is why small humps 
additionally appear in the RDF in Fig. [2], for e.g., at rn = 4.3A and at other places. 

Paniago et al. 22 have studied the stress relaxation by monitoring the RHEED pattern of 
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the growth of Pd film on Cu(lll) substrate. They found a quick increase of atom spacing 
with layer number during deposition, and the pseudomorphic growth was only up to the 
2nd film layer. Our Pd/Cu system shows consistent results. Since our substrate surface 
is the (OOl)-plane but not the (lll)-plane as used in their experiment, one more option 
exists for the film to relax the strain by transforming the structure into the close-packed 
one. We have verified other systems with smaller compressive strain, for example, the Cu/Ni 
(/ = —2.49%) system. The film layer turns to follow basically the structure of the nickel. 
These findings demonstrate the importance of lattice mismatch on thin film deposition in 
pseudomorphic growth and in determination of the structure of epilayer. 

Known that an epilayer can take a mixed structure while depositing onto a (OOl)-plane, an 
important question to be answered is to quantify the percentages of the (001)- and the (111)- 
planar structures in an epilayer and to study how they depend on the lattice mismatch. We 
propose the following method to achieve this goal. The peaks at rn = \/2a and v3a in gi{r\\) 
can serve as an unambiguous identifier to identify, respectively, the (001)- and the (111)- 
planar structures, contributed from the next nearest neighboring atoms. By integrating 
the two peaks separately with weighting 27rr\\<j£o, we calculated the mean numbers of the 
adatoms in an epilayer with a distances \/2a and \/3a, respectively, to an adatom. The 
percentages of the two structures were then computed by dividing the two mean numbers 
separately by 4 and by 6, which are, in turn, the numbers of the next nearest neighbors 
for a perfect (001)- and a perfect (lll)-structure. Subtracting these two percentages from 
1 gives the percentage for the other structure. The results, averaged from the 5th to the 
10th epilayer, which excludes the possible interface alloying region, are shown in Fig. [5] as a 
function of lattice mismatch at three temperatures 80K, 300K and 900K. 

We observed an asymmetric curve with respect to / = 0. When / is negative, the 
percentage of the (001) structure decreases enormously with increasing |/|; at the same 
moment, the percentage of the (111) structure increases. Nearly no atom sits on a site other 
than the two structures. The results shows that for a large mismatch, such as the Pt/Ni 
system (/ = —10.2%), the film layers tend to grow a close-packed structure to reduce the 
compressive strain, instead of epitaxially following the square lattice of the substrate surface. 
On the other hand, for a positive lattice mismatch, the percentage of the (OOl)-structure 
also decreases with /, but more weakly compared to the systems with negative mismatch. 
In this case, the formation of the (lll)-film structure will not be helpful to reduce the tensile 
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FIG. 5: Average percentages for the (001)-, the (111)-, and the other structures in an epilayer at 
80K, 300K and 900K as a function of lattice mismatch. 



strain; thus, almost no close-packed structure is formed. The strain is relaxed by the creation 
of breaks on the film layers. The adatoms sitting on the breaks (cf. Fig. HJ) contribute the 
percentage of the other structure. This tensile-compressive asymmetry has been emphasized 
and discussed in recent studies^. 

Concerning the effect of temperature, we found that the higher the substrate temperature, 
the less the film will grow in a pseudomorphic way; in other words, the percentage of the 
(OOl)-film structure decreases. This effect is more obvious in the regime of negative mismatch 
than positive one. For positive mismatches, the percentage of other structure also increases 
with temperature. We know that the (111) plane has the lowest surface energy for fee metals. 
Increasing temperature enhances the migration of an adatom to a more thermodynamically 
stable site, the site on the (111) plane. A negative misfit will promote this migration to form 
a close-packed structure to reduce the mechanically compressive strain, whereas a positive 
misfit will hinder it because it increases the tensile strain which is mechanically unfavorable. 
The interplay between the thermodynamics and the mechanics of the materials determines 
the film structure. 

For the implication to experiments, this study suggests that it is easier to grow pseu- 
domorphic film on (001) surface of these transition metals in the positive region of lattice 
mismatch than in the negative. We also emphasize that the method proposed in this study 
to quantify the structural composition can be adapted to experiments because the 2D RDF 
of a film layer can be experimentally obtained by surface diffraction techniques. Thus the 
structural composition versus lattice mismatch can be calculated to understand this strong 
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asymmetry predicted by our simulations. Concerning the limitation of the standard MD 
techniques, it is definitely worth in the future to investigate the effect of the neglected in- 
terlayer diffusion processes on this problem using more sophisticated simulation techniques. 

In summary, we have discussed the effect of lattice mismatch on the structure of epilayers 
on the (001) substrate surface among six transition metals. The layer coverage, the sur- 
face roughness, and the temperature effect have also been studied. We presented the first 
calculation of the composition percentage of film structure for positive and negative lattice 
mismatches. The results revealed a strong asymmetric behavior between the tensile and the 
compressive cases. The epilayer releases compressive strain by the formation of the close- 
packed planar structure, but release tensile strain by the creation of breaks and domains on 
a film layer. 
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